Z0 Data Analysis Workbook

Santiago.R, Birge Ken Tok

#import fundamental libraries
import numpy as np
import math
import scipy.constants as const
from IPython.display import IFrame
import matplotlib.pyplot as plt
import iminuit
import time
%matplotlib inline

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))

Reading Function [RD Parsons]

###Function to read in events from the ASCII data files
def read_events(fname):
    file_content=[]
    particle_content=[]
    event_content=[]
    line_content=[]
    muon_max = 0
    
    #loop over the file on a line by line basis
    for line in open(fname):
        # If the line starts with a "#" it's a comment.
        if(line.startswith("#")):
            continue
        
        # If the line has 1 entry it contains only the number of particles in the following event. 
        # Information in this line is not saved but used to define the event structures
        
        if len(line.split())==1:
            number_of_particles=int(line.split()[0])
            #initialize counter of particles and event content
            count=0
            event_content=[]
            #initialize event variables
            n_el=0
            n_mu=0
            n_had=0
            n_part=0
            passes_hadron_selection=0
            passes_muon_selection=0
            max_pt=0.
            continue
            
        #when there are more entries then it's a particle. the information about its momentum and mass are saved in a
        #structured array. another couple of variables are also added
        if len(line.split())!=1 and count<number_of_particles:
            
            # Structure containing the information about each particle within the event
            line_content=np.zeros(1, dtype = [('px','float64'),('py', 'float64'),('pz', 'float64'), # Momentum components
                                              ('m', 'float64'), # Particle mass
                                              ('p', 'float64'), # Total momentum
                                              ('e', 'float64'), # Total energy
                                              ('pt', 'float64'),('et', 'float64'), # Transverse momentum and energy
                                              ('ez', 'float64')]) # Longitudinal energy
            
            event_var=np.zeros(1, dtype = [('n_el','float64'), # Number of EM particles
                                           ('n_mu','float64'), # Number of muons
                                           ('n_had','float64'), # Number of hadrons
                                           ('n_part','float64'), # Number of particles total
                                           ('max_pt','float64'),# Maximum particle transverse momentum
                                           ('muon_max','float64') # most energetic muon REMOVE!
                                          ])
                            
            
            # Copy across the particle information from the line read in
            num, px, py, pz, m = [float(item) for item in line.split()]
            
            # Copy everything into our particle information structure
            line_content['px'] = px
            line_content['py'] = py
            line_content['pz'] = pz
            line_content['m'] = m
            line_content['p']=math.sqrt(px * px + py * py + pz * pz)
            line_content['pt']=math.sqrt(px * px + py * py)
            line_content['e']=math.sqrt(px * px + py * py + pz * pz)
            line_content['et']=math.sqrt(px * px + py * py)
            line_content['ez']=math.sqrt(pz * pz)
            event_content.append(line_content)

            if abs(abs(line_content['m'][0])-0.106)<0.001: # This is a muon
                n_mu+=1 
                if line_content['e'] > muon_max: # REMOVE
                    muon_max = line_content['e']
            elif abs(abs(line_content['m'][0])-0.000)<0.001: n_el+=1 # Photon (based on mass)
            elif abs(abs(line_content['m'][0])-0.140)<0.001: n_had+=1 # Hadron (based on mass)

            n_part=number_of_particles
            if line_content['pt'][0]>max_pt: max_pt=line_content['pt'][0]
            if count==number_of_particles-1: 
                    event_var['n_el']=n_el
                    event_var['n_mu']=n_mu
                    event_var['n_had']=n_had
                    event_var['n_part']=n_part
                    event_var['max_pt']=max_pt
                    
                    event_var['muon_max']=muon_max #REMOVE
                    muon_max=0
                    
                    file_content.append(event_var)
                    particle_content.append(event_content)

        count+=1
                    
    return np.array(file_content), particle_content

Hadronic event selection

The selected hadronic events are chosen based on two criteria;

1. The total calorimetric energy of the selected hadrons needs to be around the value of the invariant mass $\sqrt{s}$. With an energy of
   $E_{Vis}= \sum_{i=1}^n E_i = \sum_{i=1}^n \sqrt{\vec{p_i}^2+m_i^2} \approx \sum_{i=1}^n \sqrt{p_{xi}^2+p_{yi}^2+p_{zi}^2}$
   for a total of n particles, the following relation for the ratio must hold true
   $0.5 < \frac{E_{Vis}}{\sqrt{s}} < 1.5$
2. According to Chapter 5 PPE 93, the energy imbalances along the beam direction $\Delta E_{\parallel}$ and perpendicular to it $\Delta E_{\perp}$ as given by the formulas
   $\Delta E_{\perp} = \sqrt{(\sum_{i=1}^np_{xi})^2+(\sum_{i=1}^np_{yi})^2}$
   $\Delta E_{\parallel} = \sum_{i=1}^np_{zi}$
   need to also stay within a range of
   $|\Delta E_{\parallel}|/E_{vis} < 0.6$
   $\Delta E_{\perp}/E_{vis} < 0.5$

#Function for the total calorimetric energy
def vis_energy(data_particles):
    list_vis_energy = []
    for i in range(len(data_particles)):
        vis_energy = 0
        for j in range(len(data_particles[i][:])):
            vis_energy += np.sqrt(data_particles[i][j]["px"]**2+data_particles[i][j]["py"]**2+data_particles[i][j]["pz"]**2)
        list_vis_energy.append(vis_energy)
    vis_energy_array = np.array(list_vis_energy)
    return vis_energy_array

#Function for the perpendicular energy
def perp_energy(data_particles):
    list_perpendicular_energy = []
    for i in range(len(data_particles)):
        perpendicular_energy = 0
        px_n = 0
        py_n = 0
        for j in range(len(data_particles[i][:])):
            px_n += data_particles[i][j]["px"]
            py_n += data_particles[i][j]["py"]
        perpendicular_energy = np.sqrt(px_n**2+py_n**2)
        list_perpendicular_energy.append(perpendicular_energy)
    perpendicular_energy_array = np.array(list_perpendicular_energy)
    return perpendicular_energy_array

#Function for the parallel energy
def parallel_energy(data_particles):
    list_parallel_energy = []
    for i in range(len(data_particles)):
        parallel_energy = 0
        for j in range(len(data_particles[i][:])):
            parallel_energy += data_particles[i][j]["pz"]
        list_parallel_energy.append(np.abs(parallel_energy))
    parallel_energy_array = np.array(list_parallel_energy)
    return parallel_energy_array

# here we can define the cuts to make to select hadron or muon events
def make_hadron_selection(event_properties, event_particles, sqrt_s):
    #define the cuts by the above criteria
    part_cut = event_properties[:]["n_part"] >= 15       #Take Hadrons  EDIT
    criteria_cut_1_upperlim = vis_energy(event_particles)/sqrt_s < 1.5
    criteria_cut_1_lowerlim = vis_energy(event_particles)/sqrt_s > 0.5
    criteria_cut_2 = perp_energy(event_particles)/vis_energy(event_particles) < 0.5
    criteria_cut_3 = parallel_energy(event_particles)/vis_energy(event_particles) < 0.6
    # Then look for events where all cuts are passed
    selection = np.all((part_cut,criteria_cut_1_lowerlim, criteria_cut_1_upperlim,criteria_cut_2, criteria_cut_3), axis=0)
    return selection.ravel()# Ravel to make sure output is flat

hadrons_mc_properties, hadrons_mc_particles = read_events('hadrons.dat')
selection = make_hadron_selection(hadrons_mc_properties, hadrons_mc_particles , 91.2)
count = np.count_nonzero(selection)
print(count/len(hadrons_mc_properties))

0.9911

Muon Event Selection

For the selected muon events, the criteria are somewhat different than for hadronic events. First, as according to Chapter 5 PPE 93, the muon event must be within the angular range
$0<|cos(\Theta)|<0.8$ with $cos(\Theta) = \frac{p_z}{\sqrt{p_x^2+p_y^2+p_z^2}}$
Then, their energy in the form of their momentum $E_{mu} \approx \sqrt{p_{xi}^2+p_{yi}^2+p_{zi}^2}$ should be roughly $1/2\sqrt{s}$ since most of the beam energy is conserved and thus
$40GeV \leq E_{mu} \approx \sqrt{p_{xi}^2+p_{yi}^2+p_{zi}^2} \leq 50GeV$
with the total impulse of both muons lying below 5GeV since their opposite directions would make their total impulse add up to 0.
$p_{tot} = \sqrt{\vec{p_{\mu}}+\vec{p_{\overline{\mu}}}} \leq 5GeV$

#First, a function is defined that calculates the cosine of the angle theta for events within the 40-50GeV threshold. For events outside it, a value of 0 is returned,
#i.e. the output must be 0 < cos(theta) < 0.8 for the muonic event to be selected 
def cos_theta(data_particles):
    list_cos_theta = []
    for i in range(len(data_particles)):
        cos_theta = 0
        for j in range(len(data_particles[i][:])):
            if 40<= np.sqrt(data_particles[i][j]["px"]**2+data_particles[i][j]["py"]**2+data_particles[i][j]["pz"]**2) <= 50:
                cos_theta = np.abs(data_particles[i][j]["pz"]/np.sqrt(data_particles[i][j]["px"]**2+data_particles[i][j]["py"]**2+data_particles[i][j]["pz"]**2))
                break
        list_cos_theta.append(float(cos_theta))
    cos_theta_array = np.array(list_cos_theta)
    return cos_theta_array

#Next, a function is defined for the total impulse of the muon pairs
def p_tot(data_particles):
    list_p_tot = []
    for i in range(len(data_particles)):
        p_tot = np.zeros(3)
        p1 = np.zeros(3)
        p2 = np.zeros(3)
        for j in range(len(data_particles[i][:])):
            if 40<= np.sqrt(data_particles[i][j]["px"]**2+data_particles[i][j]["py"]**2+data_particles[i][j]["pz"]**2) <= 50:
                p1 = np.array([data_particles[i][j]["px"][0],data_particles[i][j]["py"][0],data_particles[i][j]["pz"][0]])
                for k in range(1,len(data_particles[i][j+1:])+1):
                    if 40<= np.sqrt(data_particles[i][k+j]["px"]**2+data_particles[i][k+j]["py"]**2+data_particles[i][k+j]["pz"]**2) <= 50:
                        p2 = np.array([data_particles[i][k+j]["px"][0],data_particles[i][k+j]["py"][0],data_particles[i][k+j]["pz"][0]])
                        break
                p_tot = p1+p2
                break
        p_tot_abs = np.sqrt(p_tot[0]**2+p_tot[1]**2+p_tot[2]**2)
        list_p_tot.append(p_tot_abs)
    p_tot_array = np.array(list_p_tot)
    return p_tot_array

def make_muon_selection(event_properties, event_particles): 
    #Particle cuts    
    part_cut = event_properties[:]["n_part"] >= 1
    part_cut_2 = event_properties[:]["n_part"] < 10
    part_cut_mu = event_properties[:]["n_mu"] >= 2
    #Muon cuts
    angle_cut_upperlim = cos_theta(event_particles) < 0.8
    angle_cut_lowerlim = cos_theta(event_particles) > 0
    impulse_cut_upperlim = p_tot(event_particles) < 10 #für mehr Muone Ereignisse kann man gerne diesen Wert auf 10GeV setzen
    impulse_cut_lowerlim = p_tot(event_particles) > 0
    selection = np.all((np.transpose(part_cut)[0], np.transpose(part_cut_2)[0], np.transpose(part_cut_mu)[0], angle_cut_upperlim, angle_cut_lowerlim, impulse_cut_upperlim, impulse_cut_lowerlim), axis=0)
    return selection.ravel()

input_data_file = "./89gev.dat" # first define which energy range data we want to read in

# Then we get our data from the input file
data_events, data_particles = read_events(input_data_file)
data_muon_selection = make_muon_selection(data_events,data_particles)
print(data_muon_selection)

[False False False ... False False False]

Select Data [RD Parsons]

input_data_file = "./89gev.dat" # first define which energy range data we want to read in

# Then we get our data from the input file
data_events, data_particles = read_events(input_data_file)

#we retrieve the muon and hadron MC data
muon_mc_events, muon_mc_particles =read_events("muons.dat")
hadron_mc_events, hadron_mc_particles =read_events("hadrons.dat")

data_var = []
print(len(data_events), "events in data, ", len(muon_mc_events), "muon MC events, ", len(hadron_mc_events),"hadron MC events")

10000 events in data,  9969 muon MC events,  10000 hadron MC events

###### MODIFIZIERT FÜR NEUE SELEKTIONSALGORYTHMEN

variable = "n_part" # The name of the variable we want to plot
xmin, xmax, nbins = 0., 50, 51 # And the binning of the histogram we want to make (min, max, number of bins)
binning_type = "linear" # And whether we want log or linear binning

# Define the figure we want to draw
figure, ((ax1, ax2, ax3)) = plt.subplots(1, 3,  figsize=(18,5.5))

# Make selections on data and store a bool of those passing
data_hadron_selection = make_hadron_selection(data_events,data_particles,89.48)
data_muon_selection = make_muon_selection(data_events,data_particles)

# Do same for MC muons
muon_mc_hadron_selection = make_hadron_selection(muon_mc_events,muon_mc_particles,91.2)
muon_mc_muon_selection = make_muon_selection(muon_mc_events,muon_mc_particles)

# and MC hadrons
hadron_mc_hadron_selection = make_hadron_selection(hadron_mc_events,hadron_mc_particles,91.2)
hadron_mc_muon_selection = make_muon_selection(hadron_mc_events,hadron_mc_particles)

# Define the bins for our histogram as an array of the bin boundaries
if binning_type is "log":
    bins = np.logspace(np.log10(xmin),np.log10(xmax), nbins)
else:
    bins = np.linspace(xmin, xmax, nbins)

# The first histogram we make is our data events vs the MC muon and hadron events with no cuts made
# When comparing histograms we weight MC so that integral is same as data
weight_muon = len(data_events) / len(muon_mc_events)
weight_hadron = len(data_events) / len(hadron_mc_events)


ax1.hist(data_events[:][variable], bins=bins, alpha=0.5, label="data")
ax1.hist(muon_mc_events[:][variable], bins=bins, alpha=0.5, label="MC muons")
ax1.hist(hadron_mc_events[:][variable], bins=bins, alpha=0.5, label="MC hadrons")
ax1.set_xscale(binning_type)
#ax1.set_yscale(binning_type)
ax1.set_ylabel("Number of events", weight="bold")
ax1.set_xlabel(variable, weight="bold")
ax1.legend()
ax1.title.set_text("No cuts")


# Next we make the same plot for events passing hadron cuts
weight_muon = np.ones(np.sum(muon_mc_hadron_selection)) * \
                (np.sum(data_hadron_selection) / np.sum(muon_mc_hadron_selection))
weight_hadron = np.ones(np.sum(hadron_mc_hadron_selection)) * \
                (np.sum(data_hadron_selection) / np.sum(hadron_mc_hadron_selection))

ax2.hist(data_events[data_hadron_selection][variable], bins=bins, alpha=0.5, label="data")
ax2.hist(muon_mc_events[muon_mc_hadron_selection][variable], weights=weight_muon,
         bins=bins, alpha=0.5, label="MC muons")
ax2.hist(hadron_mc_events[hadron_mc_hadron_selection][variable], weights=weight_hadron, 
         bins=bins, alpha=0.5, label="MC hadrons")
ax2.set_xscale(binning_type)
#ax2.set_yscale(binning_type)
ax2.set_ylabel("Number of events", weight="bold")
ax2.set_xlabel(variable, weight="bold")
ax2.legend()
ax2.title.set_text("Hadron cut")


# And finally for events passing muon cuts
weight_muon = np.ones(np.sum(muon_mc_muon_selection)) * \
                (np.sum(data_muon_selection) / np.sum(muon_mc_muon_selection))
weight_hadron = np.ones(np.sum(hadron_mc_muon_selection)) * \
                (np.sum(data_muon_selection) / np.sum(hadron_mc_muon_selection))

ax3.hist(data_events[data_muon_selection][variable], bins=bins, alpha=0.5, label="data")
ax3.hist(muon_mc_events[muon_mc_muon_selection][variable], weights=weight_muon, 
         bins=bins, alpha=0.5, label="MC muons")
ax3.hist(hadron_mc_events[hadron_mc_muon_selection][variable], weights=weight_hadron,
         bins=bins, alpha=0.5, label="MC hadrons")
ax3.set_xscale(binning_type)
ax3.legend()
ax3.title.set_text("Muon cut")

ax3.set_ylabel("Number of events", weight="bold")
ax3.set_xlabel(variable, weight="bold")

/home/santi/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:47: RuntimeWarning: divide by zero encountered in long_scalars
/home/santi/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:68: RuntimeWarning: divide by zero encountered in long_scalars

Text(0.5, 0, 'n_part')

Calculate the Hadron and Muon Cross Sections

#Calculate the cross section for hadrons

#           N_actual  - N_BKG                N_cutMC
# sigma = --------------------    epsilon = ---------
#            L   epsilon                     N_totMC

def get_epsilon(err_n):
    # Load MC data and select events
    data_MC_hadron_events, data_MC_hadron_particles = read_events("./hadrons.dat")
    data_MC_hadron_selection = make_hadron_selection(data_MC_hadron_events,data_MC_hadron_particles,91.2)
    
    data_MC_muon_events, data_MC_muon_particles = read_events("./muons.dat")
    data_MC_muon_selection = make_muon_selection(data_MC_muon_events,data_MC_muon_particles)

    count = np.count_nonzero(data_MC_hadron_selection)
    epsilon_hadron = count/len(data_MC_hadron_events)
    err_hadron = err_n[0]/len(data_MC_hadron_events)

    count = np.count_nonzero(data_MC_muon_selection)
    epsilon_mu = count/len(data_MC_muon_events)
    err_mu = err_n[1]/len(data_MC_muon_events)

    print("Epsilon_hadron =", epsilon_hadron, "+-", err_hadron)
    print("Epsilon_mu =", epsilon_mu, "+-", err_mu)
    return [epsilon_hadron,err_hadron], [epsilon_mu,err_mu]

def get_crossec_err_standard(n_prime,n_b,epsilon,L):
    sum1 = n_prime[1]/(epsilon[0]*L[0])
    sum2 = n_b[1]/(epsilon[0]*L[0])                            # err noise = 0
    sum3 = (n_prime[0]-n_b[0])*epsilon[1]/(epsilon[0]**2*L[0]) # err epsilon = 0
    sum4 = n_b[0]*L[1]/(epsilon[0]*L)
    err = np.sqrt(sum1**2+sum4**2)
    return err

def get_crossec_err(n_prime,epsilon,L): # n_b = 0
    sum1 = n_prime[1]/(epsilon[0]*L[0])
    sum2 = (n_prime[0])*epsilon[1]/(epsilon[0]**2*L[0]) # err epsilon = 0 NEEDS EDIT! GET EPSILON ERR
    sum3 = L[1]/(epsilon[0]*L[0])
    err = np.sqrt(sum1**2+sum2**2+sum3**2)
    return err

def crossec(input_data_file,integrated_luminosity,sqrt_s,est_noise,epsilons,err_n):
    # Load data and select Events
    data_events, data_particles = read_events(input_data_file)
    data_hadron_selection = make_hadron_selection(data_events,data_particles,sqrt_s)
    data_muon_selection = make_muon_selection(data_events,data_particles)

    N_hadron = np.count_nonzero(data_hadron_selection)
    N_mu = np.count_nonzero(data_muon_selection)
    print(N_hadron, N_mu)
    epsilon_hadron, epsilon_mu = epsilons

    sigma_hadron = (N_hadron-est_noise[0][0])/(integrated_luminosity*epsilon_hadron[0])
    err_sigma_hadron = get_crossec_err([N_hadron,err_n[0]],epsilon_hadron,[integrated_luminosity,integrated_luminosity*0.01])
    sigma_mu = (N_mu-est_noise[1][0])/(integrated_luminosity*epsilon_mu[0])
    err_sigma_mu = get_crossec_err([N_mu,err_n[1]],epsilon_mu,[integrated_luminosity,integrated_luminosity*0.01])
    print("Sigma_hadron =", sigma_hadron, "+-", err_sigma_hadron)
    print("Sigma_mu =", sigma_mu, "+-", err_sigma_mu)
    return [sigma_hadron,err_sigma_hadron], [sigma_mu,err_sigma_mu]
    

err_n = [25,7]
epsilon_hadron, epsilon_mu = get_epsilon(err_n)
epsilons = [epsilon_hadron, epsilon_mu]

Epsilon_hadron = 0.9911 +- 0.0025
Epsilon_mu = 0.461129501454509 +- 0.0007021767479185475

input_data_file = './89gev.dat'
integrated_luminosity = 179.3
sqrt_s = 89.48
est_noise = [[0,0] ,[0,0]] #placeholder. If modified, modify crossec() and err_crossec()
sigma_hadron_89, sigma_mu_89 = crossec(input_data_file,integrated_luminosity,sqrt_s,est_noise,epsilons,err_n)

1786 33
Sigma_hadron = 10.050407916570135 +- 0.1433048267468306
Sigma_mu = 0.39912666436232747 +- 0.08739856796811876

input_data_file = './91gev.dat'
integrated_luminosity = 135.9
sqrt_s = 91.33
est_noise = [[0,0] ,[0,0]] #placeholder. If modified, modify crossec() and err_crossec()
sigma_hadron_91, sigma_mu_91 = crossec(input_data_file,integrated_luminosity,sqrt_s,est_noise,epsilons,err_n)

3991 74
Sigma_hadron = 29.630896732204327 +- 0.20034864721187676
Sigma_mu = 1.1808353754079948 +- 0.11380045877257187

input_data_file = './93gev.dat'
integrated_luminosity = 151
sqrt_s = 93.02
est_noise = [[0,0] ,[0,0]] #placeholder. If modified, modify crossec() and err_crossec()
sigma_hadron_93, sigma_mu_93 = crossec(input_data_file,integrated_luminosity,sqrt_s,est_noise,epsilons,err_n)

2131 32
Sigma_hadron = 14.239312664168049 +- 0.17116509807975122
Sigma_mu = 0.45956836232094933 +- 0.10284534295546764

Fitting the Breit Wigner PDF [RD Parsons]

from scipy.integrate import quad
#definition of stuff you don't really need to look into, here we define the functions for the convolution between the
#breit-wigner and the gaussian
def radcorr(z,s):
    me2    = 0.0000002611200
    alphpi = 0.0023228196
    zeta2  = 1.644934
    cons1  = -2.16487
    cons2  = 9.8406
    cons3  = -6.26
    fac2   = 1.083
    l      = math.log(s/me2)
    beta   = 2.*alphpi*(l-1.)
    del1vs = alphpi*((3./2.)*l+2.*zeta2-2.)
    del2vs = pow(alphpi,2.)*(cons1*pow(l,2.)+cons2*l+cons3)
    delvs  = 1.+del1vs+del2vs
    del1h  = -alphpi*(1.+z)*(l-1.)
    delbh  = pow(alphpi,2.)*(1./2.)*pow((l-1.),2.)*((1.+z)*(3.*math.log(z)-4.*math.log(1.-z))-(4./(1.-z))*math.log(z)-5.-z)
    rad    = beta*pow((1.-z),beta-1.)*delvs+del1h+delbh
    rad *=fac2
    return rad

# Breit-Widgner distribution
def bw(s, s0, g2, m2):
    return s0*(s*g2)/(pow(s-m2,2)+m2*g2)

def bw2(x,pars):
    return bw(pow(x,2),pars[0],pow(pars[1],2),pow(pars[2],2))

def integrand2(x, p0, p1, p2, p3):
    integrand = bw(x*p0,p1,p2,p3)
    integrand *= radcorr(x,p0)
    
    return integrand

def convolution(s,s0,g2,m2):
    
    if(s<60):
        print("s too small")
        return -1

    zmin=pow(60,2)/s
    zmax=0.9999999999
    result, error = quad(integrand2, zmin,zmax, args=(s,s0,g2,m2), epsrel=1e-6)
    return result

def convolution2(x, sigma, gamma, mass):
    y_est = []
    for xp in x:
        y_est.append(convolution(pow(xp,2), sigma, pow(gamma,2), pow(mass,2)))
    return y_est

from iminuit import Minuit

#now we want to perform a fit for hadrons
#we have three energy-points
npoints=3
#here you need to declare the energy-vs-cross section and the errors from your estimation
x= [89.48,91.33,93.03]
y= [sigma_hadron_89[0],sigma_hadron_91[0],sigma_hadron_93[0]] # These values are of course example guess values, replace them with your own
ex= [0.,0.,0.]
ey= [sigma_hadron_89[1],sigma_hadron_91[1],sigma_hadron_93[1]] # These error values should also be replaced

def chi_squared(sigma, gamma, mass):
    expectation = np.array(convolution2(x, sigma, gamma, mass))
    return np.sum((expectation-np.array(y))**2/np.array(ey))

minuit = Minuit(chi_squared, sigma=70., gamma=5., mass=91.)
minuit.migrad()
minuit.hesse()
minuit.minos()

fit_parameters =  np.array([minuit.values[0], minuit.values[1], minuit.values[2]])
errors =  np.array([minuit.errors[0], minuit.errors[1], minuit.errors[2]])

# Make a plot of the best fit values in comparison to measurements
plt.figure(figsize=(6,5.5))
plt.errorbar(x, y, xerr=ex, yerr=ey, fmt="o", label="LEP Measurements")

# Evaluate our fitted function at lots of values of energy to make a line
# then draw it on top of the points
fit_x = np.linspace(80,100,100)
plt.plot(fit_x, convolution2(fit_x, *fit_parameters), color="black", linestyle="dashed", label="BW Fit")
plt.plot(fit_x, bw2(fit_x, (fit_parameters[0], fit_parameters[1], fit_parameters[2])), 
         color="red", linestyle="dotted", label="BW Fit (no correction)")

plt.xlabel("Energy (GeV)", weight="bold")
plt.ylabel("Cross section (nb)", weight="bold")
plt.legend()
plt.grid()

sigma_hadron_fit, gamma_hadron_fit, mass_hadron_fit = fit_parameters
sigma_hadron_fit_error, gamma_hadron_fit_error, mass_hadron_fit_error = errors
chisq = chi_squared(sigma_hadron_fit, np.abs(gamma_hadron_fit), mass_hadron_fit)

print(minuit.fmin)
print(minuit.params)
print("Chi^2=",chisq)

/home/santi/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:17: IMinuitWarning: errordef not set, using 1 (appropriate for least-squares)

┌──────────────────────────────────┬──────────────────────────────────────┐
│ FCN = 0.0001014                  │              Nfcn = 324              │
│ EDM = 0.000102 (Goal: 0.1)       │                                      │
├───────────────┬──────────────────┼──────────────────────────────────────┤
│ Valid Minimum │ Valid Parameters │        No Parameters at limit        │
├───────────────┴──────────────────┼──────────────────────────────────────┤
│ Below EDM threshold (goal x 10)  │           Below call limit           │
├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤
│  Covariance   │     Hesse ok     │ Accurate  │  Pos. def.  │ Not forced │
└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘
┌───┬───────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│   │ Name  │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │
├───┼───────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ sigma │   36.9    │    0.6    │    -0.6    │    0.6     │         │         │       │
│ 1 │ gamma │   2.59    │   0.07    │   -0.07    │    0.07    │         │         │       │
│ 2 │ mass  │   91.20   │   0.04    │   -0.04    │    0.04    │         │         │       │
└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
Chi^2= 0.00010143829183613398

# A simple function to convert our covariance matrix into a correlation matrix
def correlation_from_covariance(covariance):
    errors = np.sqrt(np.diag(covariance))
    outer_errors = np.outer(errors, errors)
    correlation = covariance / outer_errors
    correlation[covariance == 0] = 0
    return correlation

#correlation_matrix = correlation_from_covariance(cov_matrix)
#print(correlation_matrix)

def breit_wiegner_fit(y,ex,ey,sigma,gamma,mass):
    #now we want to perform a fit for muons, try to implement it by yourself
    #we have three energy-points
    npoints=3
    x= [89.48,91.33,93.03]
    minuit = Minuit(chi_squared, sigma, gamma, mass)
    minuit.migrad()
    minuit.hesse()
    minuit.minos()

    fit_parameters =  np.array([minuit.values[0], minuit.values[1], minuit.values[2]])
    errors =  np.array([minuit.errors[0], minuit.errors[1], minuit.errors[2]])

    # Make a plot of the best fit values in comparison to measurements
    plt.figure(figsize=(6,5.5))
    plt.errorbar(x, y, xerr=ex, yerr=ey, fmt="o", label="LEP Measurements")

    # Evaluate our fitted function at lots of values of energy to make a line
    # then draw it on top of the points
    fit_x = np.linspace(80,100,100)
    plt.plot(fit_x, convolution2(fit_x, *fit_parameters), color="black", linestyle="dashed", label="BW Fit")
    plt.plot(fit_x, bw2(fit_x, (fit_parameters[0], fit_parameters[1], fit_parameters[2])), 
         color="red", linestyle="dotted", label="BW Fit (no correction)")

    plt.xlabel("Energy (GeV)", weight="bold")
    plt.ylabel("Cross section (nb)", weight="bold")
    plt.legend()
    plt.grid()

    sigma_muon_fit, gamma_muon_fit, mass_muon_fit = fit_parameters
    sigma_muon_fit_error, gamma_muon_fit_error, mass_muon_fit_error = errors

    sigma_ret = [sigma_muon_fit,sigma_muon_fit_error]
    gamma_ret = [gamma_muon_fit,gamma_muon_fit_error]
    mass_ret = [mass_muon_fit,mass_muon_fit_error]
    chisq = chi_squared(sigma_muon_fit, gamma_muon_fit, mass_muon_fit)
    print(minuit.fmin)
    print(minuit.params)
    print("Chi^2=",chisq)
    return sigma_ret,gamma_ret,mass_ret

#we also want to see how the breit-wigner would look for muons
y= [sigma_mu_89[0],sigma_mu_91[0],sigma_mu_93[0]] # These values are of course example guess values, replace them with your own
ex= [0.,0.,0.]
ey= [sigma_mu_89[1],sigma_mu_91[1],sigma_mu_93[1]]
sigma = 70
gamma = 5
mass = 92
sigma,gamma,mass = breit_wiegner_fit(y,ex,ey,sigma,gamma,mass)

/home/santi/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:7: IMinuitWarning: errordef not set, using 1 (appropriate for least-squares)
  import sys
/home/santi/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:44: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected
  in the extrapolation table.  It is assumed that the requested tolerance
  cannot be achieved, and that the returned result (if full_output = 1) is 
  the best which can be obtained.

┌──────────────────────────────────┬──────────────────────────────────────┐
│ FCN = 2.111e-06                  │             Nfcn = 1679              │
│ EDM = 2.11e-06 (Goal: 0.1)       │                                      │
├───────────────┬──────────────────┼──────────────────────────────────────┤
│ Valid Minimum │ Valid Parameters │        No Parameters at limit        │
├───────────────┴──────────────────┼──────────────────────────────────────┤
│ Below EDM threshold (goal x 10)  │           Below call limit           │
├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤
│  Covariance   │     Hesse ok     │ Accurate  │  Pos. def.  │ Not forced │
└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘
┌───┬───────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│   │ Name  │   Value   │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit-  │ Limit+  │ Fixed │
├───┼───────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ sigma │    1.5    │    0.7    │    -0.5    │    0.7     │         │         │       │
│ 1 │ gamma │    2.3    │    1.4    │    -6.4    │    1.8     │         │         │       │
│ 2 │ mass  │   91.0    │    0.8    │    -0.8    │    0.7     │         │         │       │
└───┴───────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
Chi^2= 2.1106581644757765e-06

Computing the Partial Width of Electron $\Gamma _e$

#def muon_electron_partial_width(sigma_0,gamma_z,mass_z):
#    sigma_not_0 = sigma_0[0]*2.576e-6
#    return mass_z[0]*gamma_z[0]*np.sqrt(sigma_not_0/(12*np.pi))

def partial_width_e(T_z, M_z, sigma_0):
    part_width = np.sqrt((sigma_0*1e-6*2.576)/(12*np.pi))*T_z*M_z
    return part_width   

def muon_insecurities_partial_width(sigma,TZ,MZ, u_TZ, u_MZ, u_sigma): # I think you mean uncertainties.... have more self confidence partial width!
    u_part_sigma = MZ*TZ*np.sqrt(1/(12*np.pi*sigma))*u_sigma/2 *np.sqrt(2.576e-6)
    u_part_TZ = MZ*np.sqrt(sigma*2.576e-6/(12*np.pi))*u_TZ
    u_part_MZ = TZ*np.sqrt(sigma*2.576e-6/(12*np.pi))*u_MZ
    u_partial_width_e = np.sqrt(u_part_sigma**2+u_part_MZ**2+u_part_TZ**2)
    return u_partial_width_e

#compute the insecurities
muon_insecurities_partial_width(2.6,91.2,1.5,0.07,0.04,0.6)

0.006828673862885472

#compute the partial width of the electrons using the maximum muon cross section and lepton universality
partial_width_e(2.6,91.2,1.5)

0.07591387943682473

Computing the Decay Width of Hadrons and Neutrinos $\Gamma _{hadron}$, $\Gamma _\nu$

#compute the partial width of the neutrino
# Gam_nye =GF MZ^3/24sqrt2pi [1+(1-4|Q|sin^2theta)] 

def Gam_nye(m_z,m_zerr):
    val = 2*1.166e-5*m_z**3/(24*np.sqrt(2)*np.pi)
    err = np.sqrt((6*1.166e-5*m_z**2/(24*np.sqrt(2)*np.pi))**2*m_zerr**2)
    return [val, err]

gamma_nye = Gam_nye(mass_hadron_fit,mass_hadron_fit_error)
gamma_nye

[0.16590490611587655, 0.00020808011339244235]

# partial width of the electrons using the maximum hadron cross section and the total width
def hadron_electron_partial_width(M_z,sigma_0,Gam_Z,Gam_nye):
    a = Gam_Z[0]/3-Gam_nye[0]
    b = M_z[0]**2*sigma_0[0]*2.576e-6*Gam_Z[0]**2/(36*np.pi)
    val = 1/2*(a-np.sqrt(a**2-4*b))

    sq = np.sqrt(1/4*a**2-4*b)
    print(1/4*a,4*b,sq)
    err_Z = 1/6-(1/6*a-2*M_z[0]**2*sigma_0[0]*Gam_Z[0]/(9*np.pi))/(2*sq)
    err_N = a/(4*sq)-1/2
    err_M = M_z[0]*sigma_0[0]*Gam_Z[0]**2/(9*np.pi*sq)
    err_s = M_z[0]**2*Gam_Z[0]**2/(18*np.pi*sq)
    err=(M_z[1]/M_z[0]+sigma_0[1]/sigma_0[0]+Gam_Z[1]/Gam_Z[0]+Gam_nye[1]/Gam_nye[0])*val*4
    print(err_Z,err_N,err_M,err_s)
    print(val,err)
    #err = np.sqrt((err_Z*Gam_Z[1])**2+(err_n*Gam_nye[1])**2+(err_M*M_z[1])**2+(err_s*sigma_0[1]))
    return val,err


mass = [mass_hadron_fit,mass_hadron_fit_error]
sigma = [sigma_hadron_fit,sigma_hadron_fit_error]
gamma = [gamma_hadron_fit,gamma_hadron_fit_error]


hadron_electron_partial_width(mass,sigma,gamma,gamma_nye)

0.1742263218893954 0.18755024867763104 nan
nan nan nan nan
0.0754477302936648 0.013333478128579938

/home/santi/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in sqrt
  import sys

(0.0754477302936648, 0.013333478128579938)

Calculate Weinberg Angle $\Theta _W$

# calculate Weinberg angle
def theta_w(Gam_f,M_z,Q_f=-1):
    G_f = 1.166e-5
    a = G_f*M_z**3/(24*np.sqrt(2)*np.pi)
    print(a)
    theta_1 = (1+np.sqrt(Gam_f/a-1))/(4*np.abs(Q_f))
    theta_2 = (1-np.sqrt(Gam_f/a-1))/(4*np.abs(Q_f))
    return theta_1,theta_2

theta_w(0.0833464,91)

0.08240368155751832

(0.2767397721769793, 0.2232602278230207)

Calculate Decay Widths of Quarks $\Gamma _{u}$, $\Gamma _d$ and NC factor

G_f = 1.166*1e-5
def T_quark(Q, M_Z):
    Weinberg_Faktor = 1+(1-4*np.abs(Q)*0.22)**2
    T = G_f*(M_Z**3)/(24*np.pi*np.sqrt(2))*Weinberg_Faktor
    return T

td = T_quark(1/3, 91.2) #down quark
tu = T_quark(2/3, 91.2) #up quark

def NC(T_Hadr):
    K = 1.03
    NCol = T_Hadr/(K*(2*tu+3*td))
    return NCol

NC(1.86)

3.1829070229629925